SUBROUTINE matrix_inv3 (R1, Rinv)


! ----------------------------------------------------------------------
! SUBROUTINE: matrix_inv3.f90
! ----------------------------------------------------------------------
! Purpose:
!  3x3 Matrix inversion
! ----------------------------------------------------------------------
! Input arguments:
! - R1: 			Array 3x3
!
! Output arguments:
! - Rinv : 			Inverse matrix of input matrix 3x3
! ----------------------------------------------------------------------
! Thomas D. Papanikolaou, Geoscience Australia              30 June 2016
! ----------------------------------------------------------------------


      USE mdl_precision
      IMPLICIT NONE

! ---------------------------------------------------------------------------
! Dummy arguments declaration
! ---------------------------------------------------------------------------
! IN
      REAL (KIND = prec_q), INTENT(IN) :: R1(3,3)
! OUT
      REAL (KIND = prec_q), INTENT(OUT) :: Rinv(3,3)
! ---------------------------------------------------------------------------

! ----------------------------------------------------------------------
! Local variables declaration
! ----------------------------------------------------------------------
      REAL (KIND = prec_q) :: Det_R
      REAL (KIND = prec_q) :: Cof_R(3,3), Adj_R(3,3)
! ----------------------------------------------------------------------


! ----------------------------------------------------------------------
! Determinant of R1
      Det_R =    R1(1,1) * (R1(2,2)*R1(3,3) - R1(2,3)*R1(3,2))  & 
	           - R1(1,2) * (R1(2,1)*R1(3,3) - R1(2,3)*R1(3,1))  &
	           + R1(1,3) * (R1(2,1)*R1(3,2) - R1(2,2)*R1(3,1))
! ----------------------------------------------------------------------


! ----------------------------------------------------------------------
! Cofactor matrix
      Cof_R(1,1) =          ( R1(2,2)*R1(3,3) - R1(2,3)*R1(3,2) )
      Cof_R(1,2) = -1.0D0 * ( R1(2,1)*R1(3,3) - R1(2,3)*R1(3,1) )
      Cof_R(1,3) =          ( R1(2,1)*R1(3,2) - R1(2,2)*R1(3,1) )
	  
      Cof_R(2,1) = -1.0D0 * ( R1(1,2)*R1(3,3) - R1(1,3)*R1(3,2) )
      Cof_R(2,2) =          ( R1(1,1)*R1(3,3) - R1(1,3)*R1(3,1) )
      Cof_R(2,3) = -1.0D0 * ( R1(1,1)*R1(3,2) - R1(1,2)*R1(3,1) )

      Cof_R(3,1) =          ( R1(1,2)*R1(2,3) - R1(1,3)*R1(2,2) )
      Cof_R(3,2) = -1.0D0 * ( R1(1,1)*R1(2,3) - R1(1,3)*R1(2,1) )
      Cof_R(3,3) =          ( R1(1,1)*R1(2,2) - R1(1,2)*R1(2,1) )
! ----------------------------------------------------------------------

	  
! ----------------------------------------------------------------------
! Adjoint matrix
      Adj_R = transpose (Cof_R)
! ----------------------------------------------------------------------


! ----------------------------------------------------------------------
! Inverse matrix
      Rinv = (1.0D0 / Det_R) * Adj_R
! ----------------------------------------------------------------------



END

